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SCHIESTEL’S  DERIVATION  OF  THE  EPSILON  EQUATION  AND  TWO  EQUATION 
MODELING  OF  ROTATING  TURBULENCE 

ROBERT  RUBINSTEIN*  AND  YE  ZHOUt 


Abstract.  As  part  of  a  more  general  program  of  developing  multiple-scale  models  of  turbulence,  Schi- 
estel  suggested  a  derivation  of  the  homogeneous  part  of  the  dissipation  rate  transport  equation.  Schiestel  s 
approach  is  generalized  to  rotating  turbulence.  The  resulting  model  reproduces  the  main  features  observed 
in  decaying  rotating  turbulence. 

Key  words,  dissipation  rate  equation,  rotating  turbulence,  two-equation  turbulence  models,  multiple- 
scale  turbulence  models 


Subject  classification.  Fluid  Mechanics 

1.  Introduction.  The  dissipation  rate  transport  equation  continues  to  resist  systematic  derivation, 
either  from  the  governing  equations  or  even  from  statistical  closures.  Much  of  the  closure-based  work  is 
summarized  in  [1];  more  recent  work  is  summarized  in  [2].  In  many  respects,  the  most  successful  derivation 
of  the  e  transport  equation  is  due  to  Schiestel  [3].  Among  the  successes  of  the  derivation  is  a  rather  good 
value  Cel  =1.5  and  the  demonstration  that  necessarily,  Ce2  >  Cei> 

It  is  well-known  that  the  derivation  of  the  e  equation  in  rotating  turbulence  encounters  additional 
difficulties  because  rotation  does  not  appear  explicitly  in  the  exact  transport  equation  for  the  dissipation 
rate.  Instead,  the  effect  of  rotation  is  indirect,  entering  only  through  quantities  like  the  turbulent  time- 
scale.  In  the  present  work,  the  e  transport  equation  is  treated  by  combining  Schiestel’s  arguments  with  the 
phenomenology  for  rotating  turbulence  of  Zhou  [4].  The  most  direct  generalization  of  the  argument  of  [3] 
leads  to  a  rotation-sensitized  e  equation  with  the  same  form  as  the  standard  e  equation,  but  with  an  increased 
value  of  Ce2\  a  model  of  this  type  was  proposed  by  Okamoto  [5].  A  simple  modification  of  the  argument  of 
[3]  yields  instead  a  model  of  the  form  first  proposed  by  Bardina  et  al  [6].  The  implications  of  these  models 
for  decaying  rotating  turbulence  are  discussed. 

2.  Review  of  Schiestel’s  derivation.  We  begin  with  a  split-spectrum  model  of  high  Reynolds  number 
turbulence, 


(2.1) 


E(/c)- 


Ck^  if  k  <  kq 

a  k>  ko 


In  Eq.  (2.1),  kq  is  the  inverse  integral  scale  of  turbulence  which  marks  the  transition  between  the  inertial 
range  and  the  large  scales.  Eq.  (2.1)  is  a  special  case  of  the  models  introduced  in  [3]  in  connection  with 
multiple-scale  turbulence  models.  This  is  no  more  than  a  schematic  model  of  the  actual  energy  spectrum; 
however,  as  stressed  in  [3]  and  [7],  to  derive  a  two-equation  model,  it  is  essential  that  the  spectrum  be 


*NASA  Langley  Research  Center,  Hampton,  VA  23681-2199.  This  research  was  supported  by  the  National  Aeronautics 
and  Space  Administration  under  NASA  Contract  No.  NASl-97046  while  the  first  author  was  in  residence  at  ICASE,  NASA 
Langley  Research  Center,  Hampton,  VA  23681.  R.  Rubinstein  acknowledges  the  hospitality  and  financial  support  provided  by 
K.  Hanjalic  at  TU  Delft. 

^Lawrence  Livermore  National  Laboratory,  Livermore,  CA  94551.  This  research  was  partially  performed  under  the  auspices 
of  the  U.S.  Department  of  Energy  through  the  University  of  California  Lawrence  Livermore  National  Laboratory  under  Contract 
No.  W-7405-Eng-48. 


1 


parametrized  in  some  simple  way.  Use  of  a  more  complex  model  like  the  von  Karman  spectrum  would  lead 
to  essentially  the  same  results. 

Denote  the  energy  in  the  inertial  range  by 


and  the  energy  in  the  large  scales  by 


^0  — 


Assume  that  the  spectral  descriptors  in  Eq.  (2.1)  are  functions  of  time:  e  —  e{t)  and  kq  —  /co(^)-  It  follows 
from  Eq.  (2.2)  that 

(2.4)  k  = 

This  equation  does  not  lead  to  the  desired  e  equation  directly,  because  it  contains  the  new  unknown  kq. 

To  solve  this  problem,  we  postulate 

(2.5)  ko  = 


based  on  a  very  similar  proposal  in  [3].  In  view  of  Eq.  (2.1),  Eq.  (2.5)  is  equivalent  to 


Then  Eqs.  (2.4)-(2.6)  give 


which  can  be  re-arranged  as 


ko/Ko  = 


I  . 


with  a  rather  good  value  for  C^i  and  a  value  for  Cc2  which  depends  on  the  choice  of  p.  This  result  is 
essentially  Eq.  (27)  of  [3], 

Following  Reynolds  [8],  the  constant  p  can  be  fixed  by  appealing  to  the  behavior  of  the  large  scales  of 
motion  during  decay.  Differentiation  of  Eq.  (2.3)  gives 


and  differentiation  of  Eq.  (2.2)  gives 


Assuming  that  decay  is  self-similar,  so  that 


Eqs.  (2.9)-(2.11)  lead  as  usual  to 


^  _  3^ 

ko  Ko 


k  _2e  2  ko 

k  3  e  3/^0 


k  _  ko 
k  ko 
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corresponding,  in  Eq.  (2.8),  to  (5  =  2/9. 

It  would  seem  that  this  argument  solves  the  problem  of  deriving  the  homogeneous  e  transport  equation, 
since  it  gives  the  values  Cti  =3/2  and  Ce2  —  11/6.  But  one  can  object  that  the  assumption  Eq.  (2.6)  is 
another  way  of  stating  the  final  result:  this  equation  states  that  the  integral  scale  Kq  ^  satisfies  a  transport 
equation  in  which  the  production  term  is  absent.  Indeed,  writing 


and  substituting 


leads  to 


d 


d  3  fci/2  . 


k  =  P~e 
e  =  |[(7eiP  -  C,2e] 


which  shows  that  the  absence  of  a  production  term  in  the  length-scale  transport  equation  is  equivalent  to 
Ca  =  3/2. 

The  injection  of  energy  at  large  scales  can  certainly  cause  the  integral  scale  to  increase;  at  the  same  time, 
turbulence  production  might  be  expected  to  it  to  decrease  through  the  enhancement  of  small  scales.  Eq. 
(2.5)  states  the  dominance  of  the  first  process  over  the  second.  Although  the  validity  of  this  approximation 
is  uncertain,  the  success  of  the  argument  is  undeniable,  and  it  seems  reasonable  to  ask  what  conclusions  will 
result  if  the  same  argument  is  applied  to  another  problem. 

3.  Rotating  turbulence.  To  derive  an  e  equation  for  rotating  turbulence,  we  will  combine  the  argu¬ 
ments  of  the  previous  section  with  Zhou’s  phenomenological  model  of  rotating  turbulence  [4].  Briefly,  this 
model  postulates  that  strong  rotation  replaces  the  nonlinear  time  scale  kj e  by  the  inverse  rotation  rate  Q  ^ ; 
closure  theories  lead  to 


where  by  hypothesis,  T  ocVl  hence 


K^TE{Kf 


E{k)  =  Cgv^K-2 


For  notational  simplicity,  Q  will  denote  twice  the  absolute  value  of  the  rotation  rate  throughout. 

Adding  a  model  for  the  large  scales,  we  obtain  the  analog  of  the  split-spectrum  model  of  Eq.  (2.1)  for 
rotating  turbulence, 


E{k)  = 


Again,  we  have  the  energy  of  the  large  scales, 


Ck^  if  /c  <  /Co 

CB\/ne/c~^  if  /c  >  /Co 


^0  —  ^C'/Cq 


k  =  (7^v^/Cq  ^ 


and  the  inertial  range  energy 
(3.5) 
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Note  that  the  definition  of  the  integral  scale  implied  by  Eq.  (3.5)  differs  from  the  non-rotating  result  Eq. 

(2.2). 

Following  Schiestel,  we  differentiate  Eq.  (3.5)  to  obtain 


(3.6) 


^  ^  1  €  ko 
k  2  e  /Co 


As  before,  we  must  specify  an  equation  for  the  inverse  integral  scale  kq  in  order  to  complete  the  model. 

The  simplest  possibility  is  to  retain  Eq.  (2.5).  In  this  case,  substitution  of  the  rotation-modified  spectrum 
Eq.  (3.2)  again  leads  to  Eq.  (2.6),  but  with  a  new  constant  of  proportionality, 


(3.7)  «o/ko  “  -7e/A^ 

Following  the  previous  steps,  we  find  instead  of  Eq.  (2.8) 

(3.8)  e  =  2^P-(2  +  27)^ 

with  the  definite  prediction  that  Ca  =  2  and  Cc^  >  2. 

The  constant  7  can  be  evaluated  by  assuming  that  the  constant  in  Eq.  (2.5)  is  independent  of 
rotation.  Tentatively  accepting  the  non-rotating  result  ^  =  2/9  suggested  earlier,  Eq.  (2.5)  with  the 
rotation- dependent  energy  spectrum  Eq.  (3.2)  leads  to  7  =  2/9  and  to  the  value  C^2  =  22/9.  In  decaying 
rotating  turbulence,  Eq.  (3.8)  predicts  power-law  decay  in  time,  but  with  a  smaller  exponent  than  non¬ 
rotating  turbulence:  indeed,  following  [8],  we  have 

(3.9)  fc(t) 

and  the  increase  in  Ce2  due  to  rotation  from  11/6  to  22/9  implies  a  reduction  in  the  decay  rate. 

The  model  of  rotating  decaying  turbulence  implied  by  Eq.  (3.8)  has  been  advocated,  for  example  in 
[5],  and  more  recently  in  [9].  The  value  0^2  —  22/9  in  rotating  turbulence  can  be  compared  to  the  values 
Ce2  «  2.8  recommended  in  [5]  and  Ce2  =8/3  suggested  in  [9]. 

However,  the  available  data  is  also  consistent  with  the  conclusion  that  in  rotating  turbulence,  energy 
transfer  is  suppressed  completely,  and  energy  becomes  trapped  in  the  largest  scales  of  motion,  where  it 
undergoes  purely  viscous  decay.  This  picture,  which  is  inconsistent  with  any  kind  of  power-law  decay,  is 
advocated  for  example  by  [10]  and  [11].  Which  description  of  decaying  rotating  turbulence  is  correct  remains 
controversial;  for  now,  we  would  like  to  explore  some  models  which  are  consistent  with  the  second  viewpoint. 

The  derivation  of  Eq.  (3.2)  assumes  that  the  time-scale  in  strongly  rotating  turbulence  is  the  inverse 
rotation  rate.  This  idea  suggests  replacing  Eq.  (2.6)  by 

(3.10)  ^  =  -7'11 

in  the  strong  rotation  limit.  Eqs.  (3.6)  and  (3.10)  yield  the  e  equation  in  the  form 

(3.11)  e  =  2^{P  -  e) -'f'Cle 

The  rotation  dependence  found  in  Eq.  (3.11)  coincides  with  that  of  the  well-known  Bardina  model  [6];  we 
argued  previously  [1]  for  the  strong  rotation  limit  of  this  model  on  the  basis  of  simplified  closure  arguments. 
Integration  of  the  Bardina  model  for  decaying  turbulence  in  the  strong  rotation  limit  gives  the  results  that 
e  decays  exponentially  in  time,  but  that  the  kinetic  energy  approaches  a  constant;  if  viscosity  is  included  in 
the  analysis,  then  the  kinetic  energy  undergoes  purely  viscous  decay. 
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Although  these  conclusions  are  consistent  with  numerical  and  experimental  observations  [10],  the  as¬ 
sumption  Eq.  (3.10)  underlying  the  present  derivation  has  the  consequence  that  the  integral  scale  grows 
exponentially.  This  was  cited  in  [9]  as  evidence  against  the  Bardina  model  itself,  although  [11]  argued  that 
quite  different  two-point  behavior  can  be  consistent  with  the  same  single-point  model. 

The  difficulty  is  not  so  much  with  SchiestePs  formalism,  but  with  applying  Eq.  (3.10),  an  isotropic 
result,  to  rotating  turbulence.  In  rotating  turbulence,  the  Taylor-Proudman  theorem  forces  the  large  scales 
of  motion  to  be  nearly  two-dimensional.  Consequently,  the  integral  scales  parallel  and  perpendicular  to  the 
rotation  axis  are  unequal  [10]. 

It  is  rather  difficult  to  capture  this  effect  in  any  isotropic  model.  But  suppose  that  we  combine  Eqs. 
(3.6)  and  (2.9)  to  give 


(3.12) 

and  simply  postulate  the  large  rotation 


k  1  e  Iko 
k  2  e  3  ko 

limit  of  Eq.  (3.11)  for  decaying  turbulence 


(3.13) 


€ 

e 


Then  we  obtain 

(3.14) 


k 

k 


Ih 
3  ko 


or  equivalently, 

(3.15)  kkl^^  =  k{0)koiQy^h-'>'^*/^ 


instead  of  the  self-similarity  postulate  Eq.  (2.11)  for  non-rotating  turbulence.  Unlike  the  argument  leading  to 
Eq.  (3.8),  which  like  the  derivation  for  isotropic  turbulence  assumes  that  the  energy  decay  of  the  large  scales 
and  the  inertial  range  scales  is  linked  by  self-similarity,  the  present  derivation  instead  allows  the  dynamics 
of  the  large  scales  and  the  inertial  range  scales  to  be  different. 

The  problem  of  decaying  rotating  turbulence  is  defined  by  the  energy  equation  together  with  Eq.  (3.13) 
and  either  Eq.  (3.14)  or  Eq.  (3.15).  Numerical  integration  will  be  required  to  solve  these  equations  in 
general,  but  it  is  evident  that  these  equations  are  consistent  with  the  limits 


e  =  0 

(3.16)  k  =  0 

while 


ko  =  const. 

(3.17)  K.o  “  const. 

Thus,  the  kinetic  energy  in  the  inertial  range  vanishes,  the  energy  transfer  vanishes,  but  the  kinetic  energy 
in  the  large  scales  and  the  integral  scale  both  approach  constants  in  the  absence  of  viscosity. 

Let  us  summarize  the  differences  between  the  two  dynamic  descriptions  of  rotating  decay.  Power-law 
decay,  but  with  a  reduced  exponent,  follows  if  the  decay  of  both  the  large-scale  energy  and  the  inertial  range 
energy  is  linked  through  the  self-similarity  assumption  Eq.  (2.11).  The  alternative  description,  which  leads 
instead  to  Eqs.  (3.16)  and  (3.17)  allows  the  large-scale  and  inertial  range  energies  to  evolve  independently. 
The  argument  also  implies  that  in  the  long-time  limit,  viscous  dissipation  and  energy  transfer  are  unequal: 
energy  transfer  can  vanish,  but  viscous  dissipation  is  always  nonzero. 
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4.  Conclusions.  SchiestePs  derivation  of  the  e  transport  equation  has  been  generalized  to  rotating 
turbulence.  By  assuming  that  the  basic  scale  relationship  Eq.  (2.5)  applies  to  both  non-rotating  and  rotating 
turbulence,  we  are  led  to  the  e  equation  in  the  form  Eq,  (3.8).  This  equation  implies  algebraic  decay  in 
time  of  decaying  rotating  turbulence  with  a  smaller  decay  rate  than  non-rotating  turbulence.  Replacing 
Eq.  (3.8)  with  the  rotation-dependent  hypothesis  Eq.  (3.10)  leads  essentially  to  the  Bardina  model,  which 
implies  a  completely  different  description  of  rotating  decay:  the  nonlinear  energy  transfer  vanishes  and 
in  the  absence  of  viscous  effects,  energy  approaches  a  constant.  By  ignoring  the  two- dimensionality  and 
rotation-independence  of  the  large  scales,  this  argument  leads  to  an  incorrect  description  of  the  integral 
scale  in  decaying  rotating  turbulence.  By  modifying  SchiestePs  argument,  the  Bardina  model  is  shown  to  be 
consistent  with  saturation  of  the  integral  scale. 
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DISPLACEMENT  MODELS  FOR  THUNDER  ACTUATORS  HAVING  GENERAL  LOADS 

AND  BOUNDARY  CONDITIONS 

ROBERT  WIEMAN*,  RALPH  C.  SMITHt,  TYSON  KACKLEY^,  ZOUBEIDA  OUNAIES^  AND  JEFF  BERND^ 

Abstract.  This  paper  summarizes  techniques  for  quantifying  the  displacements  generated  in  THUN¬ 
DER  actuators  in  response  to  applied  voltages  for  a  variety  of  boundary  conditions  and  exogenous  loads. 
The  PDE  models  for  the  actuators  are  constructed  in  two  steps.  In  the  first,  previously  developed  theory 
quantifying  thermal  and  electrostatic  strains  is  employed  to  model  the  actuator  shapes  which  result  from  the 
manufacturing  process  and  subsequent  repoling.  Newtonian  principles  are  then  employed  to  develop  PDE 
models  which  quantify  displacements  in  the  actuator  due  to  voltage  inputs  to  the  piezo  ceramic  patch.  For 
this  analysis,  drive  levels  are  assumed  to  be  moderate  so  that  linear  piezoelectric  relations  can  be  employed. 
Finite  element  methods  for  discretizing  the  models  are  developed  and  the  performance  of  the  discretized 
models  are  illustrated  through  comparison  with  experimental  data. 

Key  words.  THUNDER  actuators,  shape  model,  displacement  model 

Subject  classification.  Structures  and  Materials 

1.  Introduction.  THUNDER  actuators  offer  the  capability  for  generating  large  strains  and  forces  due 
to  a  variety  of  mechanisms  including  improved  robustness  through  the  manufacturing  process  and  increased 
electromechanical  coupling  due  to  their  inherent  shape.  However,  the  full  capabilities  of  these  actuators 
have  not  yet  been  completely  quantified  either  experimentally  or  analytically  due  to  their  relatively  recent 
genesis  and  the  fact  that  their  behavior  differs  quite  substantially  from  standard  unimorphs  or  bimorphs. 
In  this  paper,  we  discuss  modeling  techniques  for  quantifying  the  displacements  generated  by  THUNDER 
actuators  in  response  to  applied  voltages  for  a  variety  of  boundary  conditions  and  exogenous  loads.  The 
development  of  corresponding  finite  element  techniques  is  also  addressed  and  the  accuracy  of  the  resulting 
finite  dimensional  models  is  illustrated  through  comparison  with  experimental  data. 

Model  development  is  considered  in  two  steps:  (i)  the  characterization  of  the  actuator  shape  as  a  function 
of  the  manufacturing  process  and  (ii)  the  development  of  a  PDE  model  for  the  actuator  behavior  based  on 
Newtonian  principles.  The  first  component  has  been  addressed  in  previous  investigations  [3]  and  only  those 
details  necessary  for  the  development  of  the  subsequent  PDE  model  will  be  discussed.  As  detailed  in  [3], 
the  characteristic  curved  shape  of  THUNDER  actuators  is  due  primarily  to  differing  thermal  coefficients 
in  the  constituent  materials,  which  produce  thermal  stresses  in  the  combined  actuator  during  cooling,  and 
secondarily  to  the  reorientation  of  dipoles  during  repoling.  We  note  that  the  quantification  of  strains  due  to 

“Department  of  Mathematics,  North  Carolina  State  University,  Raleigh,  NC  27695  (rewieman@eos.ncsu,edu).  This  research 
was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  grant  533155. 

^Department  of  Mathematics,  North  Carolina  State  University,  Raleigh,  NC  27695  (rsmith@eos.ncsu.Bdu).  This  research 
was  supported  in  part  by  the  Air  Force  Office  of  Scientific  Research  under  the  grant  AFOSR-F49620-01- 1-0107  and  in  part 
by  the  National  Aeronautics  and  Space  Administration  under  NASA  Contract  Number  NASl-97046  while  this  author  was  a 
consultant  at  ICASE,  M/S  132C,  NASA  Langley  Research  Center,  Hampton,  VA  23681-2199. 

^Department  of  Mathematics,  North  Carolina  State  University,  Raleigh,  NC  27695  (tcko.ckle@unity.ncsu.edu).  This  research 
was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  grant  533155. 

§  ICASE,  M/S  132C,  NASA  Langley  Research  Center,  Hampton,  VA  23681-2199  {z.ounaies@larc.nasa.gov).  This  research 
was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Contract  Number  NASl-97046  while  this 
author  was  in  residence  at  ICASE,  M/S  132C,  NASA  Langley  Research  Center,  Hampton,  VA  23681-2199. 

ICeramic  Materials  Laboratory  Bowen  Hall,  Room  318,  70  Prospect  Avenue,  Princeton,  NJ  08540  (JDBernd@princeton.edu). 


1 


thermal  gradients  has  been  investigated  for  a  variety  of  applications  (e.g.,  see  [2,  4,  5,  6])  with  certain  aspects 
having  been  considered  for  THUNDER  actuators  [3,  7].  The  quantification  of  stresses  and  strains  due  to 
repoling  is  based  on  domain  theory  for  general  ferroelectric  materials  [8,  9, 12,  13].  Thin  shell  theory  is  then 
employed  to  develop  PDE  models  which  quantify  the  stresses  and  displacements  throughout  the  actuator 
when  voltage  is  applied  to  the  piezoceramic  patch.  For  this  analysis,  it  is  assumed  that  the  actuators  are 
operating  at  low  to  moderate  drive  levels  for  which  linear  piezoelectric  relations  are  adequate.  Techniques 
for  extending  these  models  to  regimes  in  which  the  piezoelectric  response  is  nonlinear  and  hysteretic  are 
under  investigation  and  will  utilize  methods  outlined  in  the  concluding  remarks. 

Because  the  PDE  model  is  infinite  dimensional,  approximation  techniques  must  be  considered  to  obtain 
a  finite  dimensional  model  which  is  appropriate  for  implementation.  This  is  accomplished  through  a  hybrid 
finite  element  approach  utilizing  linear  and  cubic  Hermite  basis  functions.  This  numerical  approach  differs 
from  that  employed  in  [14],  where  a  NASTRAN  model  was  employed  to  predict  dome  heights,  in  that  the  finite 
element  method  was  developed  directly  for  the  PDE  model  used  to  characterize  the  physical  mechanisms  for 
the  actuator.  Hence  the  resulting  finite  dimensional  system  incorporates  the  physical  properties  associated 
with  the  differing  constituent  materials  thus  permitting  a  detailed  analysis  of  various  aspects  of  the  actuator 
dynamics  (e.g.,  stresses  or  strains  at  various  points  along  the  length  of  the  actuator).  This  approach  also 
permits  direct  extension  of  the  numerical  method  to  nonlinear  structural  models  for  high  drive  level  dynamics 
as  well  as  models  which  incorporate  the  hysteresis  and  constitutive  nonlinearities  inherent  to  piezoceramic 
materials  at  moderate  to  high  drive  levels. 

The  manufacturing  conditions  for  THUNDER  are  outlined  in  Section  2  and  a  model  which  quantifies 
the  resulting  curved  shape  is  summarized  in  Section  3.  The  PDE  model  quantifying  the  displacements  is 
then  presented  in  Section  4  along  with  boundary  conditions  which  characterize  a  variety  of  experimental 
setups.  The  numerical  approximation  techniques  are  discussed  in  Section  5  and  examples  illustrating  the 
performance  of  the  resulting  finite  dimensional  model  are  presented  in  Section  6.  Finally,  future  work, 
including  techniques  to  extend  the  model  to  nonlinear  and  hysteretic  regimes  will  be  outlined  in  Section  7. 

2.  Actuator  Geometry.  THUNDER  actuators  are  typically  comprised  of  a  piezoceramic  wafer,  a 
metallic  backing  material,  hot  melt  adhesive  layers,  and  an  optional  metallic  top  layer  as  depicted  in  Fig¬ 
ure  2.1a.  As  detailed  in  [3],  materials  commonly  employed  for  backing  layers  include  aluminum,  stainless 
steel  and  brass  while  LaRC-SI  is  employed  as  the  adhesive. 

During  the  manufacturing  process,  the  materials  are  placed  in  a  vacuum  bag  and  heated  to  325^^  C  under 
a  pressure  of  241.3  kPa.  During  the  cooling  process,  the  LaRC-SI  solidifies  at  approximately  270""  C  and 
subsequent  cooling  produces  curvature  in  the  actuator  due  to  differing  thermal  coefficients  of  the  constituent 
materials.  Because  the  Curie  temperature  for  PZT-5A  (350°  C)  is  in  the  proximity  of  the  manufacturing 
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Fig.  2,1.  (a)  Components  of  a  THUNDER  actuator;  (b)  Curvature  observed  in  a  THUNDER  actuator. 
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temperature,  the  final  step  in  the  fabrication  process  is  comprised  of  repoling  the  material  through  the 
application  of  a  sustained  DC  voltage. 

As  illustrated  in  Figure  2.1b,  THUNDER  actuators  have  a  characteristic  dome  shape  due  to  the  man¬ 
ifestation  of  differing  thermal  properties  in  the  PZT  and  backing  material  during  the  cooling  process  and 
the  rotation  of  dipoles  during  repoling.  In  general,  curvature  will  occur  in  both  component  directions  in  a 
rectangular  actuator;  however,  for  the  models  developed  here,  we  consider  actuators  whose  width  is  small 
as  compared  with  the  length  so  that  motion  is  predominantly  in  one  dimension.  Finally,  we  note  that  the 
curvature  in  the  actuators  occurs  only  in  regions  covered  by  the  piezoceraraic  patch  and  the  end  tabs  remain 
straight. 

3.  Model  for  the  Manufacturing  Process.  A  necessary  step  before  developing  a  model  which 
predicts  the  actuator  displacements  under  various  drive  conditions  is  the  quantification  of  the  curvature 
produced  in  the  manufacturing  process.  We  summarize  here  the  characterization  of  the  stresses  produced 
during  cooling  and  repoling  which  in  turn  produce  the  curvature.  Details  regarding  this  component  of  the 
model  can  be  found  in  [3]. 

To  accommodate  various  constructions,  we  consider  actuators  with  N  layers  and  consider  a  coordinate 
system  in  which  the  x~  z  plane  corresponds  with  the  outer  edge  of  the  backing  material  and  the  ^/-coordinate 
extends  through  the  thickness  of  the  actuator.  The  width  of  the  jth  layer  is  denoted  by  bj  while  hj  indicates 
the  thickness  of  each  layer  as  depicted  in  Figure  3.1.  The  Young’s  modulus  and  thermal  coefficient  for  the 
jth  layer  are  respectively  denoted  by  Ej  and  aj.  The  strain  at  the  outer  edge  of  the  backing  material  {y  =  0) 
is  denoted  by  eo,  and  k  denotes  the  curvature  at  the  neutral  axis.  The  change  in  temperature  during  the 
bonding  process  is  indicated  by  AT.  To  incorporate  the  strains  due  to  repoling,  it  is  also  necessary  to  employ 
the  Poisson  ratio  u  and  saturation  elect rostriction  for  PZT-5A. 

As  detailed  in  [3],  the  balancing  of  forces  and  moments  due  to  thermal  and  electrostatic  stresses  yields 
the  linear  system 


(3.1) 


AS^f 


where  S  =  [eo?  and 


(3.2) 


A  = 


/- 


Ef.i  EMhj  -  hj-i)  Ef=i  EMh]  -  hU) 

h  Ef=i  EMh]  -  hU)  -I  Ef=i  Ejbjihj  -  hU)  J 


E}=i  EjbjiajAT  -  ZI26iy\,){hj  - 


I  Ef=i  EM<^iAT  -  3/25uX,){h]  -  J 


The  Kronecker  delta,  defined  by 


1  ,  if  1/  is  in  the  piezoceramic  layer 
0  ,  otherwise 


isolates  the  electrostatic  strains  due  to  repoling  to  the  piezoceramic  layer. 

To  solve  for  eo  and  k,  and  hence  obtain  the  final  radius  of  curvature  K  —  1/k,  it  is  necessary  to 
obtain  values  for  Ej,aj  for  each  of  the  constituent  layers  in  addition  to  determining  i/Xg  for  the  PZT 
compound  being  employed.  While  the  Young’s  modulus  and  thermal  coefficients  are  catalogued  for  various 
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Fig.  3.1.  Orientation  of  the  composite  THUNDER  actuator  with  five  layers. 


PZT  compounds  and  various  metallic  backing  materials,  they  are  temperature-dependent,  and  not  easily 
quantified,  for  the  LaRC-SL  Hence  these  parameters  are  typically  estimated  through  a  least  squares  fit 
to  data  for  the  constituent  materials  and  specific  manufacturing  conditions  under  consideration.  Details 
regarding  the  validity  of  the  model  for  a  variety  of  materials  and  geometries  can  be  found  in  [3]. 

4.  Displacement  Model.  The  system  (3.1)  quantifies  the  radius  of  curvature  R  =  1/k,  for  differing 
constituent  materials,  material  dimensions,  and  manufacturing  conditions.  In  this  section,  we  develop  models 
quantifying  displacements  produced  in  the  actuators  through  the  input  of  voltages  to  the  PZT  or  applied 
loads.  Because  end  conditions  crucially  affect  the  measured  displacements,  we  consider  a  variety  of  boundary 
conditions.  Finally,  we  consider  low  drive  regimes  in  which  the  relation  between  applied  voltages  and 
generated  strains  are  approximately  linear  with  minimal  hysteresis  so  that  linear  piezoelectric  relations  can 
be  employed. 

When  modeling  the  actuator,  we  consider  two  regimes.  In  the  first,  the  entire  actuator  (including 
the  tabs)  is  assumed  to  have  the  same  initial  curvature.  This  approximation  to  the  geometry  significantly 
simplifies  the  numerical  implementation  of  the  model  but  imposes  the  assumption  that  the  tabs  are  initially 
curved.  In  the  second  configuration,  only  that  region  covered  by  PZT  is  assumed  to  be  initially  curved  and 
the  end  tabs  are  taken  to  be  flat  in  the  absence  of  an  applied  voltage  or  load.  This  accurately  represents  the 
initial  conflguration  of  the  actuator  after  the  manufacturing  process.  One  of  the  objectives  when  validating 
the  model  is  to  compare  the  performance  of  both  models  and  ascertain  tab  dimensions  when  the  latter  model 
is  sufficiently  accurate. 

4.1.  Model  1.  We  consider  initially  the  model  which  results  from  the  assumption  that  the  tabs  have 
the  same  curvature  as  that  portion  of  the  actuator  covered  by  the  PZT  patch.  The  radius  of  curvature  is 
denoted  by  R  (recall  that  R  —  1/k  can  be  predicted  using  the  model  summarized  in  Section  3)  and  the 
backing  material  is  assumed  to  have  width  &,  thickness  h,  and  Young’s  modulus  E.  The  corresponding 
parameters  for  the  PZT  layer  are  respectively  denoted  by  bpe ,  hpe  and  Epe .  The  backing  material  is  assumed 
to  extend  from  ^  =  0  to  ^  L  and  the  region  [71 , 72]  covered  by  the  patch  is  delineated  by  the  characteristic 
function 

f  1  ,  7i  <  7  <  72 

(4d)  n  • 

0  ,  otherwise 

where  71  =  R9i  and  72  =  RO2  and  ^1,^2  denote  the  angles  subtended  by  the  patch.  The  longitudinal  and 
transverse  displacements,  which  are  coupled  due  to  the  curvature,  are  respectively  denoted  by  v  and  w. 

Under  the  assumptions  of  linear  displacements,  negligible  rotational  effects  and  shear  deformation,  and 
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linear  stress-strain  relations,  force  and  moment  balancing  yields  the  static  equations 


1  dNe  . 

R  de 

1  (fMe  ,  1  *  1  (fMe 

dd^ 


Here  Nq  and  Mq  denote  the  internal  force  and  moment  resultants  and  Mq  represents  the  external  moment 
generated  by  applied  voltages  to  the  patches.  Finally,  qe  and  qn  respectively  denote  applied  longitudinal 
and  normal  loads  to  the  actuator.  As  detailed  in  [1],  the  modeling  equations  (4.2)  are  consistent  with  the 
restriction  of  modified  Donnell-Mushtari  shell  equations  to  the  actuator  geometry. 

The  internal  resultants  incorporate  the  material  properties  of  the  backing  material  and  PZT  and,  as 
derived  in  [1] ,  are  given  by 


(4.3) 


,de 


(dv  \ 


+  - 


d^w 
^~d0^ 


a2 


Xpei^) 


-Eh^  (Pw 
12R^  dO^^ 


a2 

2R 


as  d?w 

ZR?  W 


Xpe{^) 


where  aa  =  (h-/2  +  hpef'  —  (V^)^  and  as  =  (/i/2  +  hpeY  ~  (/^/2)^* 

For  low  to  moderate  drive  regimes  in  which  the  linear  piezoelectric  equations  are  sufficiently  accurate, 
the  external  moment  generated  by  the  PZT  in  response  to  an  applied  voltage  V  is 


(4.4) 


Me  =  + V)XpeW 


where  dsi  is  the  linear  piezoelectric  constant. 

We  note  that  the  model  (4.2),  with  resultants  given  by  (4.3)  and  (4.4),  neglects  material  contributions 
due  to  the  LaRC-SL  If  desired,  these  contributions  can  be  incorporated  in  the  manner  described  in  [1].  We 
also  note  that  in  the  strong  form  (4.2),  differentiation  of  the  discontinuous  material  parameters  and  patch 
inputs  yields  unbounded  components  in  the  model.  This  necessitates  the  consideration  of  an  appropriate 
weak  form  of  the  model.  As  a  prelude,  however,  it  is  necessary  to  specify  appropriate  boundary  conditions. 

We  consider  four  sets  of  boundary  conditions  which  model  the  constraints  commonly  employed  in  ex¬ 
periments:  fixed-end,  pinned-end,  sliding-end  and  free-end  conditions.  These  boundary  conditions  can  be 
applied  at  either  end  of  the  beam;  to  simplify  the  discussion,  we  summarize  them  at  the  left  end  {0  =  0)  and 
note  that  similar  expression  hold  at  0  =  L. 


(i)  Fixed-End  Conditions 


(4.5) 


v{0)  =  0 

-(0)  =  >)  =  0 


(ii)  Pinned-End  Conditions 


(4.6) 


i;(0)  —  0 

^^;(0)  =  M0{O)  =  0 
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(iii)  Sliding-End  Conditions 


it;(0)  =  i;(0)  tan((?5>i) 

(4.7)  MeiO)  0 

Ne{0)  =  --Qe{0)tM(l>c) 

(iv)  Free-End  Conditions 

(4.8)  No{0)  ==  Me{0)  =  0 

In  the  sliding  end  condition,  Qe  denotes  the  shear  force  resultant  and  (pi^^c  respectively  denote  the 
initial  angle  of  the  actuator  and  the  angle  obtained  after  a  load  is  applied  (see  Figure  4.1).  It  can  be  noted 
that  to  within  a  first-order  approximation,  (pi  and  <pc  are  related  by  the  expression  (pc  =  (Pi  Finally, 

for  implementation  purposes,  it  has  been  observed  that  physically  reasonable  results  can  be  obtained  with 
the  approximation  Q$  =  0  which  is  typically  enforced  in  first-order  shell  theory.  This  yields  the  natural 
boundary  condition  No  =  0  which  is  easily  implemented. 

We  note  that  care  must  be  exhibited  when  specifying  boundary  conditions  at  the  left  and  right  ends  of 
the  actuator  to  ensure  model  well-posedness.  For  example,  the  specification  of  free-end  conditions  at  both 
X  =  0  and  x  =  L  will  yield  rigid  body  modes  and  hence  will  not  enforce  unique  solutions  since  solutions 
differing  by  a  constant  will  be  equivalent.  For  the  experiments  reported  in  Section  6,  fixed-end  conditions 
were  enforced  at  x  -  0  and  sliding-end  conditions  were  employed  at  x  =  L. 

To  accommodate  the  discontinuities  due  to  the  patch  and  to  reduce  smoothness  requirements  on  the 
basis  functions  employed  for  numerical  approximation,  we  consider  corresponding  weak  forms  of  the  modeling 
system.  The  state  space  is  taken  to  be  X  =  x  L^(n)  where  ft  =  [0,L],  The  test  functions  depend 

upon  the  boundary  conditions  under  consideration.  For  fixed,  pinned,  or  sliding-end  boundary  conditions  at 
0  =  0  and  free-end  conditions  at  ^  =  L,  we  respectively  employ  the  spaces 

V  =  {(<p,^)  m  =  0,<p(0)  =  <p'(0)  =  0} 

(4.9)  V  =  {((P,ip)  e  X  I  m  =  0,  ^(0)  -  0} 

V  =  {{(P,<^)  e  X  \  ifio)  =  (P{0)  tM<Pi)}  • 

We  note  that  the  constraints  MeiL)  =  Ne{L)  =  0  are  natural  boundary  conditions  which  do  not  require 
any  restriction  of  the  underlying  Sobolev  spaces.  Analogous  definitions  are  employed  when  considering  other 
combinations  of  boundary  conditions. 

A  weak  form  of  the  model  is  then 

(4.10) 

for  all  {(p,  ip)  in  the  appropriate  space  V, 
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4.2.  Model  2.  A  second  model  is  derived  under  the  assumption  that  the  actuator  region  covered  by 
the  patch  is  curved  while  the  tabs  are  initially  straight.  This  model  is  constructed  by  coupling  flat  and 
curved  beam  models  through  appropriate  interface  conditions. 

To  specify  the  geometry,  we  let  7  denote  the  arclength  with  7  =  0  at  the  left  end  of  the  actuator 
as  depicted  in  Figure  4.1.  We  assume  that  the  tabs  have  equal  length  t  and  that  the  portion  covered  by 
the  piezoceramic  patch  has  arclength  5  =  (^2  “  Oi)R.  The  region  covered  by  the  patch  is  denoted  by 
flpg  =  [71,72]  t]  while  n  =  [0,L]  again  denotes  the  domain  for  the  full  actuator.  The  characteristic 

function  Xpe,  defined  in  (4.1)  delineates  that  portion  of  the  structure  covered  by  the  patch.  For  the  tabs, 
the  arclength  is  designated  by  7  —  a:  whereas  it  has  the  form  7  =  RO  in  the  curved  portion  of  the  actuator. 

For  the  curved  portion  of  the  actuator,  force  and  moment  balancing  again  yield  the  coupled  relations 

1  dNe  ^ 

R  de 

1  <fMe  ,  1  ~  1 

E?  dO^  ~  R^  de^  ' 


where  Ng,  Mg  and  Mg  are  defined  in  (4.3)  and  4.4).  The  tabs  have  infinite  radius  of  curvature  which  yields 
the  uncoupled  relations 


(4.12) 


dN^  _  ,  SM^  _  . 

da:  ’  da:2 


where  the  resultants  are  given  by 

(4.13)  N,  =  Eh~ 


Eh^  d^w 
12  dx^  ' 


Finally,  the  displacements  and  slopes  at  71  and  72  are  coupled  through  the  interface  constraints 


(4.14) 


lim  v{j)  =  iim  17(7) 

7^7^  7-»-7i‘' 


lim_  u(7)  —  lim  1^(7) 

7-+7J  7-J'7^ 


lim  w{'))  —  lim  ^^;(7) 
7 ->-71“ 


lim  u;(7)  —  lim  1^(7) 

7-^7^  7->7^ 


lim  ^(7)  = 

7->-7i 


.  dx 


dw 

dx 
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7->72 


dw 

dx 


7^7^ 


Formulation  of  the  model  in  terms  of  the  arclength  7  yields  the  unified  relations 


(4.15) 
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02  (PuJ 


Xpe(7) 
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—Eh^  dPw 
12  di^ 


02  dv  a2W 


as  cPta 


Xpeil) 


M^=  y-V(/i  +  V)xp.(7)- 

The  radius  of  curvature  is  taken  to  be  =  oo  for  the  tabs.  When  employing  the  formulation  (4.15),  it 
must  be  noted  that  the  second  derivatives  may  not  exist  at  the  points  71  and  72-  Finally,  fixed,  pinned  and 
sliding-end  boundary  conditions  are  enforced  by  employing  the  constraints  (4.5)  -  (4.7)  at  7  =  0  or  7  =  L. 

The  construction  of  the  weak  model  formulation  is  analogous  to  Model  1  with  the  exception  that  second 
derivatives  of  the  transverse  test  functions  may  not  exist  at  the  interface  points  [71,72]-  To  illustrate,  the 
space  of  test  functions  for  an  actuator  with  fixed-end  conditions  at  7  =  0  and  sliding-end  conditions  at  7  =  L 
is 

V  =  e  X  \  <t>(0)  =  0,¥5(0)  =  <^'{0)  =  0,ip{L)  =  ^{L)tani<pi)  and 

lim  lim  lim_^"(7)#  lim  tf>"(7)| . 

7->7i’  7-+72  7^72  ^ 

Analogous  spaces  are  employed  for  the  remaining  combinations  of  boundary  conditions.  A  weak  form  of  the 
model  is  then 

I  <*7  =  0 

(4.17) 

i  {""7$  +  +  *7^}  <*7  =  0 

for  all  ((^,  (p)  eV. 

5.  Numerical  Approximation  Techniques.  To  approximate  the  solution  of  (4.10)  or  (4.17),  we 
consider  Galerkin  techniques  with  basis  functions  chosen  to  satisfy  smoothness  requirements  as  well  as 
boundary  and  interface  conditions.  We  consider  first  the  system  which  arises  when  discretizing  the  model 
for  the  uniformly  curved  actuator. 


5.1.  Model  1.  We  consider  Galerkin  approximation  for  v  and  w  which  are  respectively  based  on  linear 
and  cubic  Hermite  functions.  To  define  the  bases,  consider  a  uniform  partition  of  [0,L]  with  gridpoints 
Oi  =  ih,  h  =  L/iV,  i  =  0,  •  •  • ,  A/".  For  i  =  1,  *  •  • ,  iV  -  1,  linear  splines  are  taken  to  be 


(5.1) 


{O-Oi-i) 

(Oi+i-e) 

0 


,  0  G  [0i-\ ,  0i\ 
,  otherwise  . 


The  cubic  Hermite  basis  functions  employed  to  specify  w  and  are  given  by 

r  (0-0i_i)2[2(0,-0)  +  /i] 

¥’io(^)  =  ^<!  {ei+^-0f[2{di+x-e)-h] 


0 
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e  e  [0i-\,6i] 

e  €  (6»i,6ii+i] 
otherwise 


and 

{e-ei-i)Ho-Oi)  ,  9e[ei-i,ei] 

{6i+i-$f{e-ei)  ,  ee{6i,ei+i] 

0  ,  otherwise  . 

As  detailed  in  [11],  all  three  sets  of  basis  functions  vanish  identically  outside  the  interval  The 

basis  functions  <f>o^<poOi^oi  ^-nd  (pN^^NOy^Ni  defined  similarly  on  the  intervals  [0,^i]  and 
(see  [11,  pages  49  and  57]).  The  essential  boundary  conditions  (i)-(iii)  summarized  in  (4.5)-(4.7)  are  enforced 
by  omitting  or  forming  linear  combinations  of  the  boundary  basis  functions. 

The  displacements  are  then  represented  as  linear  combinations  of  basis  functions  with  coefficients  de¬ 
termined  by  enforcing  the  constraints  provided  by  the  weak  form  of  the  model.  To  illustrate,  consider  the 
discretization  of  the  model  satisfying  fixed-end  conditions  at  ^  =  0  and  pinned-end  conditions  at  ^  =  L.  The 
approximate  solutions  are  taken  to  be 

N~1 

(5-2)  r'l 

W^{0)  =  +^'^iV>nW 

j=l  J=l 

in  the  subspace  C  V.  Analogous  expressions  are  employed  for  other  combinations  of  boundary  condi- 
tions. 

A  matrix  system  is  obtained  by  considering  the  approximate  solution  in  (4.10)  with  basis  functions 
employed  as  test  functions  (this  is  equivalent  to  projecting  the  system  (4.10)  onto  the  finite  dimensional 
subspace  H^).  This  yields  the  linear  relation 

(5.3)  Kv  =  f 


where  v  —  [t^i ,  *  *  *  >  ujv— i  ?  ? '  *  *  ?  >  *^1  ?  *  ’  *  5  ^tv]  denotes  the  vector  of  unknown  coefficients. 

5.2.  Model  2.  The  formulation  of  approximation  techniques  for  Model  2  is  accomplished  in  a  similar 
manner.  In  this  case,  we  consider  uniform  partitions  on  each  of  the  subintervals  [0,71],  [7i?72]>  [72  >^] 
define  basis  functions  on  each  subdomain.  The  displacements  5 '^3^]  ~  [^1^ 5^2  5^3  1 

are  defined  in  a  manner  analogous  to  (5.2)  with  the  interface  conditions  enforced  through  the  constraints 


=  ■y^(7i)  .  'V2  (72)  =  V3  (72) 

(71)  =  (71)  >  ‘>^2  (72)  =  W3  (72) 

dw^ ,  ,  dw^ 


dw^  /  \  dw2 
-(71)  = - 


(71) 


(72)  =  ^(72) 


dy  d7  ’  ^7  ^7 

Orthogonalization  against  the  test  functions  then  yields  a  corresponding  linear  system 


Kv  =  f 


which  can  be  solved  to  obtain  the  displacement  and  slope  coefficients  at  each  of  the  nodes  in  the  subdomains. 
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6.  Experimental  Validation.  In  this  section,  we  illustrate  the  performance  of  the  combined  model 
through  comparison  with  experimental  data.  In  the  first  example,  the  strain  model  outlined  in  Section  3  is 
used  to  predict  the  shape  of  a  variety  of  actuators  as  a  function  of  the  manufacturing  process.  This  provides 
the  initial  geometry  employed  in  Section  4  when  characterizing  displacements  generated  through  voltage 
inputs  to  the  patches.  The  performance  of  the  displacement  model  is  illustrated  in  Example  2. 

6.1.  Example  1.  Actuator  Shape 

The  model  summarized  in  Section  3  quantifies  the  thermal  and  electrostatic  strains  and  resulting  changes 
in  curvature  generated  during  the  cooling  and  repoling  of  the  actuator  during  the  manufacturing  process. 
This  provides  a  means  of  characterizing  the  radius  of  curvature  and  dome  height  for  an  actuator  as  functions 
of  properties  of  the  constituent  materials  as  well  as  the  dimensions  of  these  materials.  We  consider  in  this 
example  an  actuator  construction  comprised  of  a  stainless  steel  bottom  layer,  LaRC-SI,  PZT-5A  and  a 
protective  LaRC-SI  top  layer  (there  is  no  metallic  top  layer).  The  width  of  all  materials  was  0.5  inches  and 
the  PZT  was  1.5  inches  in  length  while  the  stainless  steel  was  2.5  inches  in  length.  Hence  5  =  1.5  and  t  —  0,5 
in  Figure  4.1.  The  PZT  was  8  mils  thick  while  the  LaRC-SI  had  a  mean  thickness  of  1  mil.  Actuators  were 
constructed  with  steel  thicknesses  ranging  from  1  mil  to  20  mils  to  illustrate  the  effect  of  backing  material 
thickness  on  the  final  dome  height  of  the  actuator.  Note  that  the  dimensions  of  the  THUNDER  devices 
considered  here  permit  the  use  of  the  1-D  model. 

The  dome  heights  h  predicted  by  the  relation  (3.1)  are  compared  with  experimental  data  from  actuators 
having  steel  thicknesses  ranging  from  3  mils  to  20  mils  in  Figure  6.1.  In  the  model,  the  parameter  values 
Ep,  =  177  X  10®  N/m^,  Esi  =  7.45  x  10®  N/m^,  Esteei  ^  173  x  10®  N/m^,  =  0.8  x  10“^  asi  =  46  x  10"® 

(23-150^  C),  asi  =  60  x  10“®  (23-150^  C),  asteei  =  9.8  x  10"®,  =  0.6  x  10"^  and  z/  =  0.3  were  estimated 

through  a  least  squares  fit  to  the  full  data  set.  It  is  observed  that  the  model  predicts  both  the  trends  and 
magnitudes  for  the  dome  heights  to  within  5%  relative  accuracy  for  all  steel  thicknesses  except  6  mils  for 
which  the  relative  error  was  8%.  This  provides  a  characterization  of  the  actuator  shape  which  can  then  be 
employed  when  modeling  subsequent  displacements  or  forces. 
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Fig.  6.1.  Model  predictions  and  experimentally  measured  dome  heights  h  as  a  function  of  steel  thickness. 


10 


6.2.  Example  2:.  Actuator  Displacements 

To  illustrate  the  performance  of  the  displacement  model  presented  in  Section  4,  we  consider  the  con¬ 
struction  described  in  Example  1  for  actuators  having  fixed-end  conditions  at  the  left  edge  (7  =  0)  and 
sliding-end  conditions  at  the  right  edge  (7  ==  L).  Model  2  was  used  to  predict  the  displacements  generated 
at  the  patch  center  =  t  +  s/2)  for  a  variety  of  steel  thicknesses  and  input  voltages. 

For  fixed  voltage  levels,  model  predictions  for  the  displacement  to  voltage  ratio  as  a  function  of  steel 
thickness  are  compared  with  experimental  data  in  Figure  6.2a.  It  is  observed  that  the  model  accurately 
predicts  the  displacement  for  actuators  having  steel  thicknesses  of  3  mils,  8  mils  and  10  mils  with  discrepancies 
observed  at  1  mil  and  6  mils.  We  note  that  the  LaRC-SI  has  the  same  thickness  as  the  steel  at  1  mil  and 
we  hypothesize  that  in  this  case,  unmodeled  viscoelastic  properties  of  the  LaRC-SI  may  be  dominating  the 
elastic  properties  of  the  steel.  The  error  in  the  6  mil  prediction  reflects  the  discrepancy  observed  in  the  dome 
height  prediction  for  that  thickness. 

To  further  illustrate  the  performance  of  the  model  at  low  drive  levels,  the  predicted  displacements  for 
input  voltage  levels  of  20  V,  80  V  and  120  V  are  compared  with  experimentally  measured  displacements  in 
Figure  6.2b.  The  steel  thickness  in  this  case  was  10  mil.  It  is  observed  that  within  this  (approximately) 
linear  range,  the  model  accurately  predicts  the  displacement  for  a  variety  of  input  levels. 

7.  Concluding  Remarks.  The  model  described  here  provides  a  technique  for  quantifying  both  the 
initial  shape  of  THUNDER  devices  due  to  the  manufacturing  process  and  displacements  generated  by  ap¬ 
plied  voltages.  The  actuator  shapes  were  modeled  through  the  quantification  of  thermal  and  electrostatic 
strains  while  Newtonian  principles  were  used  to  derive  PDE  models  characterizing  displacements  for  a  va¬ 
riety  of  boundary  conditions  and  exogenous  loads.  Both  components  of  the  model  were  illustrated  through 
comparison  with  experimental  data. 

It  should  be  noted  that  linear  theory  was  employed  when  deriving  both  components  of  the  model  and 
degradation  of  performance  is  expected  in  high  drive  regimes.  The  extension  of  both  components  to  nonlinear 
regimes  is  under  current  investigation.  The  extension  of  the  thermal  model  to  include  nonlinear  effects  is 
being  considered  in  the  context  of  theory  in  [10].  To  incorporate  the  constitutive  nonlinearities  and  hysteresis 
inherent  to  piezoceramic  materials  at  moderate  to  high  drive  levels,  the  models  developed  in  [12,  13]  are 


Fig.  6.2.  (a)  Model  predictions  and  measured  displacements  as  a  function  of  steel  thickness;  (b)  Modeled  and  measured 
displacements  as  a  function  of  input  voltage:  xxxx  (Data), -  (Model). 
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being  combined  with  the  current  model  to  accommodate  large  inputs.  Finally,  nonlinear  shell  theory  will  be 
employed  to  ascertain  limitations  in  the  linear  PDE  presented  here  when  displacements  are  large. 
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